Anthropogenic change decouples a freshwater predator’s density feedback

Intraspecific interactions within predator populations can affect predator–prey dynamics and community structure, highlighting the need to better understand how these interactions respond to anthropogenic change. To this end, we used a half-century (1969–2018) of abundance and size-at-age data from Lake Erie’s walleye (Sander vitreus) population to determine how anthropogenic alterations have influenced intraspecific interactions. Before the 1980s, the length-at-age of younger walleye (ages 1 and 2) negatively correlated with older (age 3 +) walleye abundance, signaling a ‘density feedback’ in which intraspecific competition limited growth. However, after the early 1980s this signal of intraspecific competition disappeared. This decoupling of the density feedback was related to multiple anthropogenic changes, including a larger walleye population resulting from better fisheries management, planned nutrient reductions to improve water quality and transparency, warmer water temperatures, and the proliferation of a non-native fish with novel traits (white perch, Morone americana). We argue that these changes may have reduced competitive interactions by reducing the spatial overlap between older and younger walleye and by introducing novel prey. Our findings illustrate the potential for anthropogenic change to diminish density dependent intraspecific interactions within top predator populations, which has important ramifications for predicting predator dynamics and managing natural resources.

. Predictions of how changes in the strength of intraspecific interactions could affect the relationship between predator growth (represented as body size) and population size. In (a), predator length declines as its population size increases owing to increased intraspecific interactions, such as competition or aggression, thus reducing the energy available for surplus growth. Environmental changes that (b) increase the frequency or intensity of intraspecific interactions would be expected to increase the slope of the relationship between body size and population size, with the opposite occurring if (c) the frequency or intensity of intraspecific interactions are reduced.

Results
Anthropogenic environmental change. The largest anthropogenic environmental changes occurred during the 1980s when total phosphorus (TP) inputs declined and commercial walleye harvest (also reflecting increased walleye abundances), water transparency, and temperature increased (moving from left to right primarily on NMDS axis 1 in Fig. 3a,b). Based on mean values before versus after the 1980s, TP inputs decreased from about 18,500 to 9100 metric tonnes (Fig. 2b), whereas commercial harvest increased from about 0.2 to 2.9 million kg (Fig. 2a), and water transparency increased from about 0.9 to 1.2 m ( Fig. 2c; see Supplementary Information-Sect. 2, Fig. S2.1 for changes in seasonal temperatures). We also observed increases in soluble reactive phosphorus (SRP) inputs during the 2000s, which were primarily represented by the second ordination axis (moving from bottom to top on NMDS axis 2 in Fig. 3a). Before the 2000s , mean SRP inputs from the Maumee River were about 320 metric tonnes, which then increased to about 550 metric tonnes during 2010-2018 (Fig. 2d). ; and (f) total prey abundance (CPUE; individuals·trawl min −1 ). We also plotted (a) the estimated size of the older (age 3 +) walleye population (dashed line) to illustrate its relationship to commercial harvest.

Scientific Reports
| (2023) 13:7613 | https://doi.org/10.1038/s41598-023-34408-0 www.nature.com/scientificreports/ Prey community change. The largest changes in prey-species composition began during the 1980s and further intensified after the 2000s (Fig. 3c,d). Before the 1980s (left side of NMDS axis 1 in Fig. 3c), when TP inputs were higher and commercial harvest, water transparency, and temperatures were lower, the prey community was characterized by higher abundances of prey species preferentially selected by walleye, such as emerald and spottail shiners (Notropis atherinoides and N. hudsonius), gizzard shad (Dorosoma cepedianum), and alewife (Alosa pseudoharengus; see Supplementary Information-Sect. 1, Table S1.1). After the 1980s, as TP inputs declined and commercial harvest, water transparency, and temperature increased, the prey community shifted towards dominance by less-preferred species, particularly invasive white perch (moving rightward primarily on NMDS axis 1 in Fig. 3c). This compositional change further intensified after the 2000s because composition only remained on the right side of NMDS axis 1 after this period (Fig. 3c,d). The second NMDS axis primarily represented more sporadic interannual shifts between specific species, such as years with higher abundances of freshwater drum (Aplodinotus grunniens) and gizzard shad versus higher abundances of emerald shiner and white bass (Morone chrysops; Fig. 3c).  www.nature.com/scientificreports/ Directional changes in prey-species composition were matched by coincident shifts in their trait composition, which was also primarily represented by NMDS axis 1 (Fig. 3e,f). Before the 1980s, the prey community was characterized by species that tend to reside in the water column, that are soft-rayed (i.e., no dorsal spines), that prefer cooler temperatures, and that spawn in summer (Fig. 3e). These species also tend to be planktivores or omnivores with a higher energy density. After the 1980s, community trait composition shifted towards species with an opposite set of traits. This shift intensified after the 2000s and has persisted up to 2018 (Fig. 3f). This newer prey assemblage was characterized by benthic, spiny-rayed species that prefer warmer temperatures and spawn earlier in the year (spring or spring/summer; Fig. 3e). These species also tend to be piscivores and invertivores with lower energy densities.
Changes in the walleye density feedback. The slope of the density feedback shifted from negative to positive in both age-1 and age-2 walleye (Figs. 4 and 5 and Supplementary Information-Sect. 3). This shift was related to changes in the lake environment, prey-species composition, and prey-trait composition for age-1 walleye (based on significant interactions with the respective first NMDS axes; LRTs, n = 42; Environment: L = 4.08, df = 1, R 2 lik = 0.092, P = 0.043; Species: L = 5.21, df = 1, R 2 lik = 0.12, P = 0.022; Traits: L = 7.10, df = 1, R 2 lik = 0.16, P = 0.0077; Fig. 4) and in relation to only prey-species composition for age-2 walleye (based on a significant . The marginal effects of older walleye (age 3 +) abundance on younger (age-1) walleye body size conditioned on changes in (a) the lake environment, (b) prey-species composition, and (c) prey-trait composition. Negative marginal effects indicate that the slope of the body size ~ abundance relationship is negative, specifically that younger walleye lengths tend to decline as older walleye abundance increases. A shift to a neutral then positive marginal effect indicates an inflection point where the slope changes and lengths no longer decline as abundance increases. Black lines are the estimated slope values extracted from the Generalized Least Squares models across the full range of NMDS axis values. Gray shaded areas indicate the confidence intervals of these estimates. Arrows are included in each panel at the inflection points along with which year the slopes first switched from negative to positive based on the corresponding years from the NMDS axes in Fig. 3b,d,f. www.nature.com/scientificreports/ interaction with NMDS axis 1; LRT, n = 42, L = 5.67, df = 1, R 2 lik = 0.13, P = 0.017; Fig. S3.1). We found no significant changes in relation to prey abundance, nor with any of the second NMDS axes.
Based on the conditional marginal effects, the slope of the age-1 density feedback shifted from negative to positive around a value of -0.9 on NMDS axis 1 of the environmental ordination ( Fig. 4a), around -0.3 on NMDS axis 1 of the prey-species ordination (Fig. 4b), and around -0.5 on NMDS axis 1 of the prey-trait ordination (Fig. 4c). These values corresponded to changes during the late 1970s to mid-1980s in aspects of the lake environment (specifically TP, harvest, temperature, and water transparency) and prey (e.g., losses of alewife and shiners versus gains in white perch), which were primarily associated with NMDS axis 1 in their respective ordinations. Similarly, the slope of the age-2 density feedback also shifted from negative to positive around a value of -0.6 on NMDS axis 1 of the prey-species ordination (Fig. S2), which corresponded to shifts in prey-species composition during the late 1970s to early 1980s (Fig. 3d).
Our breakpoint analyses further supported the early 1980s as the timing of major changes in the density feedback. Based on the mid-1980s as the likely timing of the shift in the age-1 and age-2 density feedbacks, we selected 1980 through 1990 as the candidate years for our breakpoint analyses. These models then indicated that 1982 was the best-supported breakpoint year for the shift in the density feedback for both age-1 and age-2 walleye (Fig. 5).

Discussion
Our long-term (1969-2018) analyses of the Lake Erie ecosystem indicated that a shift occurred in the slope of the walleye density feedback, from negative to positive, likely during the early 1980s. The nature of this shift indicated a decoupling of the dependence of younger walleye lengths on older walleye abundance, signaling a reduction in the intensity of intraspecific competition. This shift was related to human-driven changes in commercial walleye harvest, nutrient inputs (primarily TP), water transparency, and temperature, in addition to a shift in both preyspecies and prey-trait composition. These findings raise two key questions: (1) Why did the signal of intraspecific competition disappear? and (2) What mechanisms generated the recent positive relationship? Figure 5. Relationships between the mean total length of (a) age-1 and (b) age-2 walleye and the estimated total number of older (age 3 +) walleye in Lake Erie during 1974-2015. Years progress across a color gradient from yellow to blue, with orange and purple representing intermediate decades. Solid lines indicate the best pre-breakpoint relationships, whereas dashed lines represent post-breakpoint relationships (modeled using piecewise linear regression). Gray lines show relationships with a downweighted influence of potentially high leverage points (age- 1: 1974, 1983, 1984, and 2007; age-2: 1974, 1983, 1984, and 2008). Relationships were initially negative (age-1: slope = − 51.2, R 2 lik = 0.37, P = 0.041; age-2: slope = − 42.0, R 2 lik = 0.36, P = 0.046), but switched to positive after the breakpoint (age-1: slope = 52.2, R 2 lik = 0.13, P = 0.028; age-2: slope = 46.4, R 2 lik = 0.12, P = 0.042). Note the log-scale on the x-axis. www.nature.com/scientificreports/ One explanation for the disappearance of the density feedback in the early 1980s is that anthropogenic environmental change reduced competition by creating greater spatial separation between younger and older walleye. All age-classes of walleye generally congregate in the western basin, and thus interact, during the spring spawning season 22 . After spawning, many older walleye migrate to the central and eastern basins or other connected lakes (e.g., Lake Huron) where water temperatures are cooler and the abundance of preferred prey is higher 22 . The extent of these migrations may have increased owing to multiple abiotic and biotic drivers. One potential driver is population size. Walleye migrations into other basins and nearby lakes, such as Lake Huron, can be more extensive during high abundance years (e.g., 2008-2009 35 ) versus low abundance years (e.g., 2011-2014 36 ). The dramatic anthropogenic increase in the size of the walleye population during the 1980s (also reflected by an increase in commercial harvest) may therefore have driven more extensive migrations of older walleye out of the western basin. Lake Erie has also been warming since at least the 1960s 37 and a warmer western basin may be causing older walleye to leave earlier in the season 38 . This inference is supported by reports of earlier migrations of older walleye in recent decades, with migrations mostly finishing by June during the 2010s 39 compared to July during the 1990s 40 , suggesting that migrations may have been occurring progressively earlier as the lake warmed. Additionally, declines in TP inputs during the 1980s may have furthered these migrations by improving hypolimnetic dissolved oxygen conditions and thus habitat availability in the central and eastern basins 41,42 . However, reduced nutrient inputs have likely played a smaller role compared to population size and temperature given that nutrient inputs and its impairments have increased in Lake Erie since the mid-1990s 43 without an associated change in the density feedback. Ultimately, although we hypothesize that the shift in the density feedback could be driven by changing walleye migrations, we cannot confirm this possibility due to a lack of historical walleye movement data.
Changes in the prey community may also have contributed to walleye movements and thus a reduction in intraspecific competition. Older walleye migrate not only to seek out better thermal habitat but also in pursuit of preferred prey 39,40 . The replacement of preferred prey in the western basin (i.e., pelagic, soft-rayed, and higher energy species such as shiners and alewife 44 ) with less-preferred species (i.e., benthic, spiny-rayed, lower energy invaders like white perch) could therefore contribute to earlier or more extensive migrations of older walleye in pursuit of better prey opportunities. This shift in prey composition was likely driven by several anthropogenic changes to the ecosystem. First, declines in nutrient inputs owing to phosphorus abatement programs 41 have been associated with reduced zooplankton biomass 45 and thus losses of zooplanktivorous prey preferred by walleye 27,44 . Second, better fishery and nutrient management have respectively led to increased walleye abundance 22 and water transparency 41 , which contributed to declines in preferred prey by increasing predation pressure 44,46 . Finally, invasive white perch expanded dramatically during the early 1980s owing to reduced competition (via declines in other prey species 47 ) and climate warming that reduced the mortality of white perch during cold winters 48 . Collectively, this combination of planned and unplanned anthropogenic changes likely drove the shift from preferred to less-preferred prey in Lake Erie.
An alternative (or additional) explanation for the shift in the density feedback is that the increase in white perch provided a new prey resource for walleye, thus reducing intraspecific competition. The timing of the initial shift in the density feedback around the early 1980s coincides closely with the expansion of invasive white perch in Lake Erie 49 . Our dataset first records this species in 1977 and its abundance increased by an order of magnitude in 1980 and again in 1984 and it has dominated the western basin prey community ever since (on average comprising over 50% of total prey abundance annually). Other notable Lake Erie invaders, such as round goby and dreissenid mussels, did not proliferate until after the early 1990s. It may seem counterintuitive for white perch to reduce intraspecific competition for prey given that total prey abundance has not changed (see Fig. 2f) and white perch is both well-defended by its dorsal spines and has one of the lowest energy densities of the prey community (Supplementary Information -Sect. 4). However, even well-defended, low energy prey can be beneficial for walleye, which is an opportunistic generalist that can rapidly adjust to shifts in prey composition 46,50 and can feed and grow well on prey communities dominated by less-preferred species 51 . Additionally, white perch may provide a novel resource by occupying a different temporal and spatial niche compared to walleye's preferred prey, which could augment overall prey availability 33,52 . For instance, walleye and white perch both spawn during spring, whereas walleye's preferred prey typically spawn during summer. The early life stages of white perch may therefore fill a temporal gap in prey availability for walleye. In support of this assertion, walleye diets in Lake Erie can include large quantities of white perch larvae 27 and walleye in other lakes prey heavily on white perch during spring until preferred prey become more abundant later in the year 53 . White perch also occupy habitats lower in the water column (benthopelagic or benthic) compared to the primarily pelagic preferred prey. Invasions of other benthic fishes have increased prey availability for walleye (e.g., round goby 34,50 ) and we suspect the same could be true for white perch. This novel benthic resource may therefore have reduced intraspecific competition for pelagic prey, thus contributing to the decoupling of the density feedback.
We cannot fully explain the emergence of the positive relationship between older walleye abundance and younger walleye lengths in recent decades, even after conducting several post-hoc analyses (detailed below). For example, we explored the possibility that the relationship is neutral and only appears positive because it is driven by a few extreme values in certain years, specifically 1974, 1983, 1984, 2007, and 2008. However, downweighting the influence of these potentially high leverage years did not alter our results (Fig. 5). We also investigated the possibility that intraspecific competition now affects aspects of younger walleye other than length, such as weight or maturation. While we only have data on these factors starting in the 1980s, we again found no indication of a negative relationship between these metrics and older walleye abundance (Supplementary Information -Sect. 5), providing further evidence that intraspecific competition truly has relaxed (although note that competition within age classes is still important 54 ). Finally, if length is no longer influenced by intraspecific competition between younger and older individuals, then a positive relationship could result if both metrics now respond in the same way to another component of the ecosystem. Prey availability is one possibility because high prey  26 . We found support for this explanation in a shift to a more positive relationship between younger walleye lengths and prey abundance after 1982 (Supplementary Information -Sect. 6, Fig. S6.1), which matched the estimated timing of the shift in the density feedback. This finding suggests that younger walleye lengths may now be more related to prey availability rather than intraspecific competition with older individuals. Additionally, environmental factors may be exerting stronger control compared to earlier decades now that the influence of intraspecific competition has weakened. For instance, walleye lengths during recent decades are related to lake water levels 54 , which can determine the availability and quality of littoral spawning habitat 55 . Because high-quality spawning conditions can also improve walleye recruitment 56 , years with suitable water levels could translate to both larger young walleye and higher abundances of older individuals. Our results have important implications for both basic and applied science. From a basic science perspective, our results provide evidence that density feedbacks within a species may be more plastic than expected 14,15 and can shift in ecologically surprising ways in response to anthropogenic environmental change (i.e., changing unpredictably or counterintuitively 57 ). We consider our results surprising because Lake Erie's walleye population exhibited a signal of intraspecific competition between older and younger walleye during the 1970s, whereas during the 1980s when the population became an order of magnitude larger, competition should have intensified but instead it diminished. Walleye also now counterintuitively seem to be larger during years with more individuals in the lake. While we can only speculate on the exact mechanism(s) underlying this finding, it provides novel evidence that expected intraspecific interactions can weaken in response to anthropogenic change and the nature of established relationships can even be reversed. Similar ecological surprises are emerging from novel terrestrial and aquatic ecosystems worldwide 57 , suggesting that our example may just be one of many counterintuitive changes that have occurred in predator populations. Such changes must be identified and studied to predict future changes in predator population dynamics, predator-prey interactions, and associated communities.
From an applied science perspective, our results also show the need for resource managers to periodically revisit assumptions made regarding density feedbacks in managed populations. Harvest allocation strategies for maintaining wildlife populations at sustainable levels often involve population models that incorporate density feedbacks 11,58 . Such models frequently assume that increasing harvest in years of high population abundance will benefit population growth by reducing intraspecific interactions with younger individuals. For example, in Lake Erie, allowing for high walleye harvest in years of high population abundance has been suggested to be beneficial to walleye growth and recruitment to the fishery because it improves resource availability for younger walleye 27 . However, harvesting more older individuals may now have little effect on the growth of younger individuals given that the density feedback no longer exists. Sudden shifts in the dynamics of key species can also be indicative of broader population-or ecosystem-level changes that may need to be incorporated into management efforts 59,60 . Continuing with Lake Erie as an example, the disappearance of the density feedback is itself indicative of broader changes in the environment and prey community, including potential changes in walleye movements and their impacts on the food web (e.g., feeding elsewhere and supplementing their diet with an invasive species). These examples demonstrate how management strategies may need to be altered, or new strategies developed, to account for anthropogenic changes in expected density feedbacks and the impacts of these changes on the broader ecosystem. Only by monitoring and periodically re-evaluating the relationship between predator growth and predator population size can we develop appropriate strategies to prevent detrimental outcomes, sustainably manage exploited populations, and conserve threatened ones.

Methods
We focused our analyses on Lake Erie's western basin, which is where younger and older walleye generally interact during spring when the older individuals aggregate for spawning 22 . During summer and fall, older walleye tend to migrate eastwards into the lake's deeper central and eastern basins in pursuit of cooler temperatures and more abundant prey, whereas the younger individuals tend to remain in the western basin 22 . Lake Erie environment. We collected data on four aspects of the lake environment that have been impacted by humans and that can affect walleye and their intraspecific interactions: commercial harvest, nutrient availability, water transparency, and temperature.
Commercial harvest -We quantified walleye harvest levels using lakewide data on the annual total commercial walleye harvest (kg) from US and Canadian waters during 1969-2018 from the Great Lakes Fishery Commission 61 . Most of this harvest has occurred on the Canadian side since the closure of the US commercial walleye fishery in 1970 22 . Commercial harvest reflects both the number of walleye removed from Lake Erie and (to a degree) long-term anthropogenic trends in walleye population size given that harvest tended to be low during the 1970s when the walleye population was smaller and higher during the 1980s and onwards as the population recovered (Fig. 2a).
Nutrient availability and water transparency -We quantified nutrient availability using annual inputs of total phosphorus (TP; metric tonnes) during 1969-2018 62 (sourced from www. bluea ccoun ting. org). To capture changes in more bioavailable forms of phosphorus, we also collected data on annual inputs of soluble reactive phosphorus (SRP; metric tonnes) from the Maumee River during 1975-1978 and 1982-2018 (source: www. ncwqr. org/ monit oring). Although SRP data were not available for the unmonitored years (specifically 1969-1974 and 1979-1981), SRP is positively related to river discharge 63 and river discharge data were available. Therefore, we estimated SRP for the unmonitored years based on a linear regression of SRP inputs to river discharge during 1975-1995 (R 2 = 0.44; P = 0.002).
Because walleye is a visual predator and prey use low light levels as a predation refuge 32 , we also included an index of water transparency (Secchi depth to the nearest 10 cm) to account for nutrient-driven changes in the light environment 41 Supplementary Information -Sect. 1). Fall Secchi depth is the most complete water transparency dataset available for Lake Erie and reflects long-term trends in the light environment caused by changes in phosphorus inputs 41 . Temperature -We represented changes in Lake Erie temperatures using a combination of: (i) modeled monthly mean surface water temperatures during 1969-2018, obtained from the National Oceanic and Atmospheric Administration (NOAA) Great Lakes Seasonal Hydrological Forecasting System; and (ii) daily mean minimum observed air temperatures during 1969-2018 measured at the Toledo, OH airport (www. ncdc. noaa. gov), which were highly correlated to maximum temperatures and were included as observed data to support modeled water temperatures. We converted monthly and daily values to annual seasonal means by averaging temperatures for winter (December-February), spring (March-May), summer (June-August), and fall (September-November) in each year and then converted these variables to z-scores by centering them to their respective means and dividing by their standard deviations. We then used Principal Component Analysis (PCA) to decompose the scaled seasonal surface water and air temperatures into their principal axes of variation and extracted the first PCA axis for our analyses because it represents most of the annual variability in temperature (59.8%; Supplementary Information -Sect. 2).
Prey community. To quantify prey abundance and composition, we used annual catch-per-unit-effort (CPUE; individuals·min -1 ) data from daytime bottom trawls conducted by the ODNR-DOW in the western basin in fall (mid-September through mid-October) during 1969-2018. We calculated annual prey-species composition based on the abundances of 11 species common in walleye diets (Supplementary Information -Sect. 1, Table S1.1). We also collated data on prey functional traits because intraspecific interactions among walleye may respond more strongly to changes in prey traits rather than species. We selected traits that might respond to environmental changes in Lake Erie or that are relevant to interactions with walleye. These traits included physiological tolerances to anthropogenic stress, habitat and feeding preferences, defense mechanisms, and prey energy density. All trait values were obtained from existing databases and the primary literature (Supplementary Information -Sect. 4).
Walleye body size and abundance. Annual length-at-age data for younger walleye (age-1 and age-2) were available during 1974-2015 from ODNR-DOW gillnet surveys conducted during fall typically after thermal stratification has ended and when walleye have accomplished most of their growth for the year 27 . The number of individuals caught in these surveys ranged from 148 to 4268 depending on the size of the spawning walleye population, with an annual average of about 1300 ± 100 individuals (mean ± SD). However, no age-1 walleye were caught in 2003 owing to low reproduction so we predicted the mean length for this year based on a positive linear relationship to the mean length of the age-2 population in the following year (R 2 = 0.64; P < 0.001).
The abundances of older (age 3 +) walleye during 1969-2018 were based on lakewide annual estimates of the total size of the age 3 + population obtained from statistical catch-at-age model outputs produced by the Lake Erie Committee's Walleye Task Group 64 . Statistical catch-at-age models are a type of integrated population model, which incorporate multiple data sources, including fishery-dependent and fishery-independent data, into a single analysis to estimate population demographics and dynamics 65 . Two different statistical catch-at-age models have been used to produce the estimates of the size of the adult walleye population: an initial catch-atage model used during 1969-1999 and a modern version using data from 1978 to present 64,66 . To create a single time-series for 1969-2018, we combined the 1969-1977 estimates with the 1978-2018 estimates and further corrected the initial catch-at-age values based on a linear regression of the overlapping years between both models (1978-1999; R 2 = 0.95, P < 0.001). Note that no live fishes were used in this study (only data) and all fish collections by the ODNR-DOW were conducted in adherence to the American Fisheries Society's "Guidelines for the Use of Fishes in Research" (https:// fishe ries. org/ policy-media/ scien ce-guide lines/ guide lines-for-the-useof-fishes-in-resea rch/).

Statistical analyses.
To quantify anthropogenic changes in the lake, we decomposed the environmental variables, prey taxonomic composition, and prey trait composition into their principal axes of variation using Non-metric Multidimensional Scaling (NMDS) with two-dimensional solutions performed via the 'vegan' package in R 67 . Before this analysis, the environmental variables were converted to z-scores, prey-species composition was converted to Hellinger-transformed abundances, and prey-trait composition was converted to community-weighted mean trait values 68 .
To test for changes in the density feedback, we used eight Generalized Least Squares (GLS) regression models (performed via the 'nlme' package in R 69 ) to determine whether the slope of the relationship between younger walleye body size and older walleye abundance changed in relation to changes in the lake environment or prey community. We used separate models because age-1 and age-2 walleye may respond differently to competition with older walleye 70 , and because we found that changes in the lake environment and prey community cooccurred and were thus correlated (see Fig. 3) so they could not be examined in the same models. All models related either the mean annual length of age-1 (four models) or age-2 (four models) walleye to the annual abundance of older walleye, which was also log-transformed to improve linearity. To test for changes in the slope of this relationship, we then included interactions with the annual values for either the lake environment (two NMDS axes), total prey abundance (summed CPUE across all species; Fig. 2f), prey-species composition (two NMDS axes), or prey-trait composition (two NMDS axes). Annual values do not fully represent the conditions that younger walleye have experienced across their two (age-1) or three (age-2) years of life, thus we also averaged all predictors across each year and the previous year for the age-1 models and each year and the previous www.nature.com/scientificreports/ two years for the age-2 models. Each model also included a predictor of the annual number of growing degree days during the gill net survey to control for interannual differences in survey duration and a first-order autoregressive structure to control for temporal autocorrelation (Supplementary Information -Sect. 5). All predictors were converted to z-scores prior to modeling. Significant (P < 0.05) interactions were determined by removing them from the model and comparing the fuller versus reduced models using Likelihood Ratio Tests (LRTs) and maximum likelihood. We also report the partial R 2 of these relationships (R 2 lik ; calculated using the 'R2.lik' function from the 'rr2' package in R 71 ), which is the difference in the variance explained between the model with and without the interaction based on maximum likelihood.
To visualize any detected interactions, we calculated the marginal effect of older walleye abundance on younger walleye body size conditioned on the effect of any interacting predictor (performed via the 'marginaleffects' package in R 72 ). The result is a series of estimated slopes and their confidence intervals for the body size ~ abundance relationships extracted from each GLS model across the full range of values of any interacting NMDS axes. We also determined the approximate timing of any major shifts in these relationships based on which NMDS axis drove the interaction and the timing of change along this axis that corresponded to an inflection in the slope of the density feedback (e.g., shifting from negative to neutral). In summary, the GLS models tested whether the slope of the body size ~ abundance relationship changed in relation to the lake environment or prey community and the marginal effects analyses visualize the nature of these slope changes.
While the GLS and marginal effects analyses determined how the density feedback changed, we needed an additional analysis to verify the estimated timing of these changes. To do so, we conducted iterative piecewise regressions 73 (also called segmented regression or 'breakpoint' analysis) of the relationships between older walleye abundance and age-1 and age-2 walleye body size. This method is used to identify where an inflection occurs in the relationship between x and y, with the iterative part referring to iteratively testing a set of candidate breakpoints to identify which breakpoint produces the best-fitting models. We used this analysis to fit models to earlier versus later data given that we found the body size ~ abundance relationship changed through time. We chose potential breakpoints of 1980-1990 based on the marginal effects analyses that suggested the relationship changed around the early 1980s. We fit new GLS models to the data from before and after these different candidate years and selected the breakpoint with the lowest mean squared error across models. Significance (P < 0.05) of breakpoint models was also determined using LRTs, with R 2 lik for these models reflecting the variance explained by the term for older walleye abundance.